source('../env.R')
It seems reasonable to expect that cities with simialr regional pools will have similar species entering the city, and thus a similar response to urbanisation.
communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))
Rows: 2462 Columns: 9── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, jetz_species_name, seasonal, presence, origin
dbl (1): city_id
lgl (3): present_urban_high, present_urban_med, present_urban_low
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
communities_summary = communities %>% group_by(city_id) %>% summarise(
regional_pool_size = n(),
urban_size_high = sum(present_urban_high),
urban_size_med = sum(present_urban_med),
urban_size_low = sum(present_urban_low)
)
communities_summary_long = bind_rows(
communities_summary %>% rename(urban_pool_size = 'urban_size_high') %>% dplyr::select(city_id, regional_pool_size, urban_pool_size) %>% mutate(is_urban_threshold = 'High'),
communities_summary %>% rename(urban_pool_size = 'urban_size_med') %>% dplyr::select(city_id, regional_pool_size, urban_pool_size) %>% mutate(is_urban_threshold = 'Medium'),
communities_summary %>% rename(urban_pool_size = 'urban_size_low') %>% dplyr::select(city_id, regional_pool_size, urban_pool_size) %>% mutate(is_urban_threshold = 'Low')
)
communities_summary_long
city_points = st_centroid(read_sf(filename(CITY_DATA_OUTPUT_DIR, 'city_selection.shp')))
Warning: st_centroid assumes attributes are constant over geometries of xWarning: st_centroid does not give correct centroids for longitude/latitude data
community_data_metrics = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'community_assembly_metrics.csv')) %>%
dplyr::select(city_id, mntd_normalised, fd_normalised, urban_pool_size, is_urban_threshold) %>%
left_join(read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv'))) %>%
left_join(communities_summary_long) %>%
left_join(city_points[,c('city_id', 'city_nm')]) %>%
rename(city_name='city_nm') %>%
mutate(is_urban_threshold = factor(is_urban_threshold, levels = c('low', 'medium', 'high'), labels = c('Low', 'Medium', 'High'), ordered = T)) %>%
arrange(city_id)
Rows: 792 Columns: 93── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): is_urban_threshold
dbl (92): city_id, mntd_normalised, mntd_actual, mntd_min, mntd_max, mntd_mean, mntd_sd, fd_normalised, fd_actual, fd_min, fd_max, fd_mean, fd_sd, FRic_normal...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rows: 342 Columns: 2── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id, urban_pool_size, is_urban_threshold)`Joining with `by = join_by(city_id)`
community_data_metrics
Load trait data
traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_jetz.csv'))
Rows: 304 Columns: 5── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): jetz_species_name
dbl (4): gape_width, trophic_trait, locomotory_trait, mass
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(traits)
fetch_normalised_traits = function(required_species_list) {
required_traits = traits %>% filter(jetz_species_name %in% required_species_list)
required_traits$gape_width_normalised = normalise(required_traits$gape_width, min(required_traits$gape_width), max(required_traits$gape_width))
required_traits$trophic_trait_normalised = normalise(required_traits$trophic_trait, min(required_traits$trophic_trait), max(required_traits$trophic_trait))
required_traits$locomotory_trait_normalised = normalise(required_traits$locomotory_trait, min(required_traits$locomotory_trait), max(required_traits$locomotory_trait))
required_traits$mass_normalised = normalise(required_traits$mass, min(required_traits$mass), max(required_traits$mass))
traits_normalised_long = required_traits %>% pivot_longer(cols = c('gape_width_normalised', 'trophic_trait_normalised', 'locomotory_trait_normalised', 'mass_normalised'), names_to = 'trait', values_to = 'normalised_value') %>% dplyr::select(jetz_species_name, trait, normalised_value)
traits_normalised_long$trait = factor(traits_normalised_long$trait, levels = c('gape_width_normalised', 'trophic_trait_normalised', 'locomotory_trait_normalised', 'mass_normalised'), labels = c('Gape Width', 'Trophic Trait', 'Locomotory Trait', 'Mass'))
traits_normalised_long
}
fetch_normalised_traits(c('Aplopelia_larvata', 'Chalcophaps_indica', 'Caloenas_nicobarica'))
Read in our phylogeny
phylo_tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
ggtree(phylo_tree, layout='circular')
Load resolve ecoregions
resolve = read_resolve()
Warning: st_buffer does not correctly buffer longitude/latitude datadist is assumed to be in decimal degrees (arc_degrees).
Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees
to_species_matrix = function(filtered_communities) {
filtered_communities %>%
dplyr::select(city_id, jetz_species_name) %>%
distinct() %>%
mutate(present = TRUE) %>%
pivot_wider(
names_from = jetz_species_name,
values_from = "present",
values_fill = list(present = F)
) %>%
tibble::column_to_rownames(var='city_id')
}
community_nmds = function(filtered_communities) {
species_matrix = to_species_matrix(filtered_communities)
nmds <- metaMDS(species_matrix, k=2, trymax = 30)
nmds_result = data.frame(scores(nmds))
nmds_result$city_id = as.double(rownames(scores(nmds)))
rownames(nmds_result) = NULL
nmds_result
}
https://www.datacamp.com/tutorial/k-means-clustering-r
scree_plot = function(community_nmds_data) {
# Decide how many clusters to look at
n_clusters <- min(10, nrow(community_nmds_data) - 1)
# Initialize total within sum of squares error: wss
wss <- numeric(n_clusters)
set.seed(123)
# Look over 1 to n possible clusters
for (i in 1:n_clusters) {
# Fit the model: km.out
km.out <- kmeans(community_nmds_data[,c('NMDS1','NMDS2')], centers = i, nstart = 20)
# Save the within cluster sum of squares
wss[i] <- km.out$tot.withinss
}
# Produce a scree plot
wss_df <- tibble(clusters = 1:n_clusters, wss = wss)
scree_plot <- ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
geom_point(size = 4) +
geom_line() +
geom_hline(linetype="dashed", color = "orange", yintercept = wss) +
scale_x_continuous(breaks = c(2, 4, 6, 8, 10)) +
xlab('Number of clusters')
scree_plot
}
cluster_cities = function(city_nmds, cities_community_data, centers) {
set.seed(123)
kmeans_clusters <- kmeans(city_nmds[,c('NMDS1', 'NMDS2')], centers = centers, nstart = 20)
city_nmds$cluster = kmeans_clusters$cluster
cities_community_data %>% left_join(city_nmds) %>% mutate(cluster = as.factor(cluster))
}
plot_nmds_clusters = function(cluster_cities) {
cluster_cities %>% dplyr::select(city_id, city_name, NMDS1, NMDS2, cluster) %>% distinct() %>%
ggplot(aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point() + geom_label_repel(aes(label = city_name))
}
plot_city_cluster = function(city_cluster_data_metrics, title) {
species_in_cluster = communities %>%
filter(city_id %in% city_cluster_data_metrics$city_id) %>%
group_by(jetz_species_name, city_name, present_urban_high, present_urban_med, present_urban_low) %>%
summarise() %>%
mutate(present_score = present_urban_high + present_urban_med + present_urban_low)
species_in_cluster$present_score = factor(species_in_cluster$present_score, levels = c(0, 1, 2, 3), labels = c('Regional Pool', 'Low Threshold', 'Medium Threshold', 'High Threshold'))
tree_cropped <- ladderize(drop.tip(phylo_tree, setdiff(phylo_tree$tip.label, species_in_cluster$jetz_species_name)))
gg_tree = ggtree(tree_cropped)
gg_presence = ggplot(species_in_cluster, aes(x=city_name, y=jetz_species_name)) +
geom_tile(aes(fill=present_score)) + scale_fill_manual(values = c("green", "yellow", "orange", "red")) +
theme_minimal() + xlab(NULL) + ylab(NULL) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
labs(fill='Presence')
species_in_cluster_traits = fetch_normalised_traits(species_in_cluster$jetz_species_name)
gg_traits = ggplot(species_in_cluster_traits, aes(x = trait, y = jetz_species_name, size = normalised_value)) + geom_point() + theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text.y=element_blank()) + xlab(NULL) + ylab(NULL) + labs(size = "Normalised Value")
gg_cities_mntd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = mntd_normalised, fill = is_urban_threshold)) + geom_bar(stat = "identity", position = position_dodge(preserve = "single")) + scale_fill_manual(values = c("yellow", "orange", "red")) + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("MNTD") + ylim(0, 1)
gg_cities_fd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = fd_normalised, fill = is_urban_threshold)) + geom_bar(stat = "identity", position = position_dodge(preserve = "single")) + scale_fill_manual(values = c("yellow", "orange", "red")) + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("FD") + ylim(0, 1)
gg_title = ggplot() + labs(title = title) + theme_minimal()
gg_presence %>% insert_top(gg_cities_mntd, height = 0.5) %>% insert_top(gg_cities_fd, height = 0.5) %>% insert_left(gg_tree, width = 0.75) %>% insert_right(gg_traits, width = 0.5) %>% insert_top(gg_title, height = 0.06)
}
REGION_DEEP_DIVE_FIGURES_OUTPUT = mkdir(FIGURES_OUTPUT_DIR, 'appendix_regional_deep_dive')
nearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Nearctic')
nearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "San Jose" "Los Angeles" "Concord" "Tijuana" "Bakersfield"
[6] "Fresno" "Sacramento" "Mexicali" "Hermosillo" "Las Vegas"
[11] "Phoenix" "Tucson" "Durango" "Portland" "Chihuahua"
[16] "Aguascalientes" "Seattle" "Ciudad Juárez" "San Luis PotosÃ" "Mexico City"
[21] "Saltillo" "Vancouver" "Albuquerque" "Monterrey" "Nuevo Laredo"
[26] "San Antonio" "Austin" "Houston" "Dallas" "Oklahoma City"
[31] "Calgary" "New Orleans" "Kansas City" "Omaha" "St. Louis"
[36] "Bradenton" "Tampa" "Minneapolis [Saint Paul]" "Atlanta" "Orlando"
[41] "Louisville" "Chicago" "Indianapolis" "Milwaukee" "Virginia Beach"
[46] "New York"
nearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% nearctic_cities_community_data$city_id))
Run 0 stress 0.1161938
Run 1 stress 0.1152814
... New best solution
... Procrustes: rmse 0.03018951 max resid 0.1758497
Run 2 stress 0.1152814
... New best solution
... Procrustes: rmse 0.00004112701 max resid 0.0002340076
... Similar to previous best
Run 3 stress 0.2542273
Run 4 stress 0.1448163
Run 5 stress 0.1198631
Run 6 stress 0.1497522
Run 7 stress 0.162113
Run 8 stress 0.1152814
... New best solution
... Procrustes: rmse 0.000008203033 max resid 0.00002624385
... Similar to previous best
Run 9 stress 0.1145839
... New best solution
... Procrustes: rmse 0.04424235 max resid 0.1777742
Run 10 stress 0.1602924
Run 11 stress 0.1148409
... Procrustes: rmse 0.04440444 max resid 0.1744703
Run 12 stress 0.1180524
Run 13 stress 0.1176404
Run 14 stress 0.1148408
... Procrustes: rmse 0.04441218 max resid 0.1745631
Run 15 stress 0.1145839
... New best solution
... Procrustes: rmse 0.000007210488 max resid 0.00001868172
... Similar to previous best
Run 16 stress 0.1180524
Run 17 stress 0.1161938
Run 18 stress 0.1148409
... Procrustes: rmse 0.04440316 max resid 0.1744719
Run 19 stress 0.114161
... New best solution
... Procrustes: rmse 0.005102649 max resid 0.02607172
Run 20 stress 0.1152814
Run 21 stress 0.114161
... Procrustes: rmse 0.000005376852 max resid 0.00001889389
... Similar to previous best
*** Solution reached
nearctic_cities_nmds
scree_plot(nearctic_cities_nmds)
nearctic_cities = cluster_cities(city_nmds = nearctic_cities_nmds, cities_community_data = nearctic_cities_community_data, centers = 5)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(nearctic_cities)
nearctic_biomes = st_crop(resolve[resolve$REALM == 'Nearctic',c('REALM')], xmin = -220, ymin = 0, xmax = 0, ymax = 70)
although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries
ggplot() +
geom_sf(data = nearctic_biomes, aes(geometry = geometry)) +
geom_sf(data = nearctic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_clusters.jpg'))
Saving 7.29 x 4.51 in image
nearctic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Neartic cluster 1')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster1.jpg'))
Saving 7.29 x 4.51 in image
nearctic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Neartic cluster 2')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster2.jpg'))
Saving 7.29 x 4.51 in image
nearctic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Neartic cluster 3')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster3.jpg'))
Saving 7.29 x 4.51 in image
nearctic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Neartic cluster 4')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster4.jpg'))
Saving 7.29 x 4.51 in image
nearctic_cities %>% filter(cluster == 5) %>% plot_city_cluster('Neartic cluster 5')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster5.jpg'))
Saving 7.29 x 4.51 in image
neotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Neotropic')
neotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Culiacán" "Guadalajara" "Morelia" "Acapulco" "Querétaro"
[6] "Cuernavaca" "Puebla" "Oaxaca" "Xalapa" "Veracruz"
[11] "Tuxtla Gutiérrez" "Villahermosa" "Guatemala City" "San Salvador" "San Pedro Sula"
[16] "Mérida" "Tegucigalpa" "Managua" "San José" "Cancún"
[21] "Guayaquil" "Chiclayo" "Panama City" "Trujillo" "Quito"
[26] "Havana" "Cali" "Lima" "Pereira" "Miami"
[31] "MedellÃn" "Ibagué" "Cartagena" "Kingston" "Bogota"
[36] "Barranquilla" "Bucaramanga" "Cúcuta" "Maracaibo" "Arequipa"
[41] "Barquisimeto" "Santo Domingo" "Maracay" "El Alto [La Paz]" "Caracas"
[46] "Cochabamba" "Viña del Mar [ValparaÃso]" "RÃo Piedras [San Juan]" "Barcelona" "Concepción"
[51] "Santiago" "Mendoza" "Salta" "Cordoba" "Asuncion"
[56] "Buenos Aires" "La Plata" "Ciudad del Este" "Montevideo" "Mar del Plata"
[61] "Porto Alegre" "São Paulo" "Santos" "Sao Jose dos Campos"
neotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% neotropic_cities_community_data$city_id))
Run 0 stress 0.134619
Run 1 stress 0.141222
Run 2 stress 0.1348238
... Procrustes: rmse 0.01072062 max resid 0.05476751
Run 3 stress 0.1348237
... Procrustes: rmse 0.01071747 max resid 0.05475478
Run 4 stress 0.134433
... New best solution
... Procrustes: rmse 0.006824883 max resid 0.04575085
Run 5 stress 0.134619
... Procrustes: rmse 0.006831176 max resid 0.04573137
Run 6 stress 0.1405634
Run 7 stress 0.1348046
... Procrustes: rmse 0.006552356 max resid 0.04609573
Run 8 stress 0.141222
Run 9 stress 0.1346226
... Procrustes: rmse 0.007200092 max resid 0.04572692
Run 10 stress 0.1414222
Run 11 stress 0.1363758
Run 12 stress 0.1363675
Run 13 stress 0.1344331
... Procrustes: rmse 0.00002617473 max resid 0.00009028585
... Similar to previous best
Run 14 stress 0.1344509
... Procrustes: rmse 0.00118731 max resid 0.007067449
... Similar to previous best
Run 15 stress 0.134619
... Procrustes: rmse 0.006829458 max resid 0.04574559
Run 16 stress 0.134433
... Procrustes: rmse 0.00002418546 max resid 0.00008384028
... Similar to previous best
Run 17 stress 0.134433
... New best solution
... Procrustes: rmse 0.000008862921 max resid 0.00003722846
... Similar to previous best
Run 18 stress 0.141222
Run 19 stress 0.134433
... Procrustes: rmse 0.00003797808 max resid 0.0001268688
... Similar to previous best
Run 20 stress 0.1348046
... Procrustes: rmse 0.006553689 max resid 0.04610243
*** Solution reached
neotropic_cities_nmds
scree_plot(neotropic_cities_nmds)
neotropic_cities = cluster_cities(city_nmds = neotropic_cities_nmds, cities_community_data = neotropic_cities_community_data, centers = 5)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(neotropic_cities)
neotropic_biomes = resolve[resolve$REALM == 'Neotropic',c('REALM')]
ggplot() +
geom_sf(data = neotropic_biomes, aes(geometry = geometry)) +
geom_sf(data = neotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_clusters.jpg'))
Saving 7.29 x 4.51 in image
neotropic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Neotropic cluster 1')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster1.jpg'))
Saving 7.29 x 4.51 in image
neotropic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Neotropic cluster 2')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster2.jpg'))
Saving 7.29 x 4.51 in image
neotropic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Neotropic cluster 3')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster3.jpg'))
Saving 7.29 x 4.51 in image
neotropic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Neotropic cluster 4')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster4.jpg'))
Saving 7.29 x 4.51 in image
neotropic_cities %>% filter(cluster == 5) %>% plot_city_cluster('Neotropic cluster 5')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster5.jpg'))
Saving 7.29 x 4.51 in image
palearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Palearctic')
palearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Lisbon" "Porto" "Marrakesh" "Seville" "Dublin" "Málaga"
[7] "Madrid" "Glasgow" "Bilbao" "Liverpool" "Bristol" "Manchester"
[13] "Birmingham" "Leeds" "Newcastle upon Tyne" "Sheffield" "Nottingham" "Valencia"
[19] "London" "Toulouse" "Paris" "Barcelona" "Rotterdam [The Hague]" "Brussels"
[25] "Amsterdam" "Lyon" "Marseille" "Dusseldorf" "Nice" "Frankfurt am Main"
[31] "Zurich" "Oslo" "Stuttgart" "Hamburg" "Genoa" "Nuremberg"
[37] "Copenhagen" "Munich" "Berlin" "Dresden" "Rome" "Prague"
[43] "Stockholm" "Poznan" "Vienna" "Wroclaw" "Zagreb" "Gdansk"
[49] "Budapest" "Krakow" "Warsaw" "Helsinki" "Riga" "Belgrade"
[55] "Lviv" "Sofia" "Thessaloniki" "Minsk" "Athens" "Kyiv"
[61] "Istanbul" "Odesa" "Samsun" "Luxor" "Tel Aviv" "Jerusalem"
[67] "Tbilisi" "Yerevan" "Abu Dhabi" "Dubai" "Bishkek"
palearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% palearctic_cities_community_data$city_id))
Run 0 stress 0.04203818
Run 1 stress 0.05064867
Run 2 stress 0.07206189
Run 3 stress 0.04240445
... Procrustes: rmse 0.03948386 max resid 0.2064496
Run 4 stress 0.1164521
Run 5 stress 0.05197845
Run 6 stress 0.1332585
Run 7 stress 0.05064862
Run 8 stress 0.08350482
Run 9 stress 0.04984194
Run 10 stress 0.065365
Run 11 stress 0.07340577
Run 12 stress 0.04203776
... New best solution
... Procrustes: rmse 0.00009496175 max resid 0.0003265871
... Similar to previous best
Run 13 stress 0.05819636
Run 14 stress 0.05530089
Run 15 stress 0.1107012
Run 16 stress 0.04223506
... Procrustes: rmse 0.04136181 max resid 0.2169081
Run 17 stress 0.1026199
Run 18 stress 0.04223585
... Procrustes: rmse 0.04135834 max resid 0.2172361
Run 19 stress 0.08634193
Run 20 stress 0.07085973
*** Solution reached
palearctic_cities_nmds
scree_plot(palearctic_cities_nmds)
Warning: did not converge in 10 iterationsWarning: did not converge in 10 iterationsWarning: did not converge in 10 iterationsWarning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3550)Warning: did not converge in 10 iterationsWarning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3550)
palearctic_cities = cluster_cities(city_nmds = palearctic_cities_nmds, cities_community_data = palearctic_cities_community_data, centers = 6)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(palearctic_cities)
palearctic_biomes = resolve[resolve$REALM == 'Palearctic',c('REALM')]
ggplot() +
geom_sf(data = palearctic_biomes, aes(geometry = geometry)) +
geom_sf(data = palearctic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_clusters.jpg'))
Saving 7.29 x 4.51 in image
palearctic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Palearctic cluster 1')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster1.jpg'))
Saving 7.29 x 4.51 in image
palearctic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Palearctic cluster 2')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster2.jpg'))
Saving 7.29 x 4.51 in image
palearctic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Palearctic cluster 3')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster3.jpg'))
Saving 14 x 6 in image
palearctic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Palearctic cluster 4')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster4.jpg'))
Saving 7.29 x 4.51 in image
palearctic_cities %>% filter(cluster == 5) %>% plot_city_cluster('Palearctic cluster 5')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster5.jpg'))
Saving 7.29 x 4.51 in image
palearctic_cities %>% filter(cluster == 6) %>% plot_city_cluster('Palearctic cluster 6')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster6.jpg'))
Saving 7.29 x 4.51 in image
afrotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Afrotropic')
afrotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Cape Town" "Johannesburg" "Pretoria" "Kigali" "Kampala" "Arusha" "Nairobi" "Addis Ababa"
afrotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% afrotropic_cities_community_data$city_id))
Run 0 stress 0.06095916
Run 1 stress 0.06095916
... New best solution
... Procrustes: rmse 0.00001059092 max resid 0.00001521966
... Similar to previous best
Run 2 stress 0.1254289
Run 3 stress 0.126785
Run 4 stress 0.126785
Run 5 stress 0.06095935
... Procrustes: rmse 0.0001624829 max resid 0.0003126969
... Similar to previous best
Run 6 stress 0.3007162
Run 7 stress 0.126785
Run 8 stress 0.126785
Run 9 stress 0.06095923
... Procrustes: rmse 0.0006753449 max resid 0.001228519
... Similar to previous best
Run 10 stress 0.06095916
... Procrustes: rmse 0.000006448791 max resid 0.00001242034
... Similar to previous best
Run 11 stress 0.06095923
... Procrustes: rmse 0.00004902177 max resid 0.00008420442
... Similar to previous best
Run 12 stress 0.1246142
Run 13 stress 0.126785
Run 14 stress 0.126785
Run 15 stress 0.1973021
Run 16 stress 0.1246142
Run 17 stress 0.06095904
... New best solution
... Procrustes: rmse 0.0002457047 max resid 0.0004426671
... Similar to previous best
Run 18 stress 0.0609592
... Procrustes: rmse 0.0002974389 max resid 0.0005449775
... Similar to previous best
Run 19 stress 0.1354466
Run 20 stress 0.126785
*** Solution reached
afrotropic_cities_nmds
scree_plot(afrotropic_cities_nmds)
afrotropic_cities = cluster_cities(city_nmds = afrotropic_cities_nmds, cities_community_data = afrotropic_cities_community_data, centers = 4)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(afrotropic_cities)
afrotropic_biomes = resolve[resolve$REALM == 'Afrotropic',c('REALM')]
ggplot() +
geom_sf(data = afrotropic_biomes, aes(geometry = geometry)) +
geom_sf(data = afrotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_clusters.jpg'))
Saving 7.29 x 4.51 in image
afrotropic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Afrotropic cluster 1')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster1.jpg'))
Saving 7.29 x 4.51 in image
afrotropic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Afrotropic cluster 2')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster2.jpg'))
Saving 7.29 x 4.51 in image
afrotropic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Afrotropic cluster 3')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster3.jpg'))
Saving 7.29 x 4.51 in image
afrotropic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Afrotropic cluster 4')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster4.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities_community_data = community_data_metrics %>% filter(core_realm == 'Indomalayan')
indomalayan_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Swat City" "Srinagar" "Jamnagar" "Jammu" "Rajkot" "Bikaner" "Jodhpur"
[8] "Jalandhar" "Ahmedabad" "Bhavnagar" "Ludhiana" "Anand" "Udaipur" "Surat"
[15] "Vadodara" "Ajmer" "Chandigarh" "Vasai-Virar" "Mumbai" "Jaipur" "Delhi [New Delhi]"
[22] "Nashik" "Dehradun" "Kota" "Pune" "Haridwar" "Dhule" "Ujjain"
[29] "Indore" "Ahmadnagar" "Kolhapur" "Jalgaon" "Agra" "Aurangabad" "Sangli"
[36] "Belagavi" "Gwalior" "Budaun" "Bareilly" "Dharwad" "Bhopal" "Bhind"
[43] "Mangaluru" "Solapur" "Vijayapura" "Akola" "Latur" "Kannur" "Davanagere"
[50] "Thalassery" "Amravati" "Kalaburagi" "Kozhikode" "Guruvayur" "Malappuram" "Lucknow"
[57] "Thrissur" "Mysuru" "Kochi" "Nagpur" "Kollam" "Jabalpur" "Ettumanoor"
[64] "Hyderabad" "Coimbatore" "Bengaluru" "Thiruvananthapuram" "Tiruppur" "Faizabad" "Erode"
[71] "Prayagraj" "Pratapgarh" "Salem" "Dindigul" "Madurai" "Tiruchirappalli" "Durg"
[78] "Vellore" "Tirupati" "Raipur" "Bilaspur" "Vijayawada" "Puducherry" "Chennai"
[85] "Kathmandu" "Colombo" "Rajamahendravaram" "Patna" "Kandy" "Bihar Sharif" "Visakhapatnam"
[92] "Ranchi" "Brahmapur" "Jamshedpur" "Darjeeling" "Siliguri" "Cuttack" "Bhubaneshwar"
[99] "Jalpaiguri" "Berhampore" "Kolkata" "Krishnanagar" "Guwahati [Dispur]" "Agartala" "Silchar"
[106] "Dimapur" "Bangkok" "George Town" "Kuala Lumpur" "Phnom Penh" "Singapore" "Hong Kong"
[113] "Sha Tin" "Hsinchu" "Taichung" "New Taipei [Taipei]" "Tainan" "Denpasar" "Kaohsiung"
[120] "Kota Kinabalu"
indomalayan_cities_nmds = community_nmds(communities %>% filter(city_id %in% indomalayan_cities_community_data$city_id))
Run 0 stress 0.1199339
Run 1 stress 0.1221765
Run 2 stress 0.1574272
Run 3 stress 0.165272
Run 4 stress 0.1180044
... New best solution
... Procrustes: rmse 0.03780687 max resid 0.2489361
Run 5 stress 0.1223128
Run 6 stress 0.128065
Run 7 stress 0.1462967
Run 8 stress 0.1177488
... New best solution
... Procrustes: rmse 0.01878365 max resid 0.1301728
Run 9 stress 0.1480994
Run 10 stress 0.1205201
Run 11 stress 0.1380584
Run 12 stress 0.1377785
Run 13 stress 0.1241019
Run 14 stress 0.1174802
... New best solution
... Procrustes: rmse 0.003927371 max resid 0.03013146
Run 15 stress 0.1617252
Run 16 stress 0.1645131
Run 17 stress 0.1603814
Run 18 stress 0.1241049
Run 19 stress 0.1596718
Run 20 stress 0.1228515
Run 21 stress 0.1597143
Run 22 stress 0.1547365
Run 23 stress 0.1488411
Run 24 stress 0.1464657
Run 25 stress 0.1441819
Run 26 stress 0.1221274
Run 27 stress 0.1195283
Run 28 stress 0.115795
... New best solution
... Procrustes: rmse 0.02711525 max resid 0.2451229
Run 29 stress 0.1224208
Run 30 stress 0.1167315
*** No convergence -- monoMDS stopping criteria:
2: no. of iterations >= maxit
27: stress ratio > sratmax
1: scale factor of the gradient < sfgrmin
indomalayan_cities_nmds
scree_plot(indomalayan_cities_nmds)
indomalayan_cities = cluster_cities(city_nmds = indomalayan_cities_nmds, cities_community_data = indomalayan_cities_community_data, centers = 6)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(indomalayan_cities)
indomalayan_biomes = resolve[resolve$REALM == 'Indomalayan',c('REALM')]
ggplot() +
geom_sf(data = indomalayan_biomes, aes(geometry = geometry)) +
geom_sf(data = indomalayan_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_clusters.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities %>% filter(cluster == 1) %>% plot_city_cluster('Indomalayan cluster 1')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster1.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities %>% filter(cluster == 2) %>% plot_city_cluster('Indomalayan cluster 2')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster2.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities %>% filter(cluster == 3) %>% plot_city_cluster('Indomalayan cluster 3')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster3.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities %>% filter(cluster == 4) %>% plot_city_cluster('Indomalayan cluster 4')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster4.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities %>% filter(cluster == 5) %>% plot_city_cluster('Indomalayan cluster 5')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster5.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities %>% filter(cluster == 6) %>% plot_city_cluster('Indomalayan cluster 6')
`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster6.jpg'))
Saving 7.29 x 4.51 in image